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We consider the critical properties of points of continuous glass transition as one can find in liquids 
in presence of constraints or in liquids in porous media. Through a one loop analysis we show that 
the critical Replica Field Theory describing these points can be mapped in the (^''-Random Field 
Ising Model. We confirm our analysis studying the finite size scaling of the p-spin model defined on 
sparse random graph, where a fraction of variables is frozen such that the phase transition is of a 
continuous kind. 

pv^ PACS numbers: 05.20.-y, 75.10Nr 

o 

CN I. INTRODUCTION 

o 

^^ The last years of research have emphasized the importance of fluctuations in understanding glassy phenomena. The 

^H present comprehension of long lived dynamical heterogeneities in supercooled liquids compares the growth of their 

^^ typical size to the appearance of long range correlations at second order phase transition points [T] . Unfortunately, in 

^Sl supercooled liquids, the theoretical study of these correlations beyond the Mean Field is just at an embryonic level. 

It has been recently proposed that the putative discontinuous dynamical transition of Mode Coupling Theory, which 

is present when all activated processes are neglected, belongs to the universality class of the unstable cj)^ theory in a 

Q random field ((/)^-RFIM in the following) [2_. However in real systems the activated processes cannot be neglected, 

1^ the only remnant of the transition is a dynamical crossover and it is not clear if there is a range where the prediction 

G of the theory can be tested. 

I In usual phase transition often we have a line of first order transitions that ends at a second order terminal critical 

jrt point. The most popular case are ferromagnets: at low temperatures there is a first order transition when the magnetic 
field crosses zero (the magnetization has a discontinuity) and this transition lines end at the usual critical point. The 
same phenomenon happen for the gas liquid transition: it is a first order transition at low temperatures that ends in 
C^ a second order transition at the critical point. 

S A similar situation can occur for liquids undergoing a glass transition, where lines of discontinuous glass transitions 

I can terminate in critical points. In this note we focus our attention to these terminal points, where the glass transition 

■ O becomes continuous and activation does not play a major role in establishing equilibrium. This transition has both a 

^ dynamic and a thermodynamic character, and it is not necessarily wiped out in finite dimension. Glassy critical points 

7^ have been theoretically studied in detail both at the dynamic and at the thermodynamic level. In dynamical Mode 

I— I Coupling Theory (MCT)[3j these are known as A^ singularities, and have been recently observed in simulations of 

kinetically constrained models on Bethe lattice [1] . At the thermodynamic level they are known from mean field Spin 

\^ Glass models [S] and Integral Equations approximations of liquid theory [B] . At these points the discontinuity in the 

JiT Edwards- Anderson non-ergodicity parameter vanishes and correspondingly, the separation of dynamics in alpha and 

^b beta regime is blurred. On approaching the critical points from the discontinuous transition side the MCT exponents 

QQ characterizing the beta relaxation go to zero and the alpha relaxation follows a universal scaling function j7]. In 

■^ general bulk liquids display a discontinuous transition pattern. However, the transition can become continuous for a 

-^ particular choice of the parameters, e.g. in presence of constraints or of quenched disorder. It has been argued that 

f~^ within MCT the glass transition can become continuous for liquids are confined in porous media 8J . If one studies 

^Sj the transition as a function of the spatial density of the confining matrix pM , one finds lines of discontinuous dynamic 

T-H and thermodynamic transition that get displaced at lower and lower liquid densities, until they merge at a common 

critical point where the transition becomes continuous. 

From the theoretical side a suitable way of constraining a glassy system consists in introducing a "pinning field" term 
in the Hamiltonian pushing the system in the direction of a randomly chosen reference equilibrium configuration. In [5] 
H it was proposed a phase diagram in the plane of temperature and pinning field, showing lines of first order dynamical 

and thermodynamical transition that merge and terminate in a common critical point as reproduced in figure 1. This 
view, based on simple spin glass models, was confirmed for liquids in the replica hypernetted chain approximation in 
[TO] and supported by numerical simulation of realistic model liquids in (TUl HJ . 

More recently, the interest for the phase diagram of constrained systems has been renewed by liquid simulations 
where a finite fraction of the particles are frozen to the position they take in a selected equilibrium configuration[l2] . 
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FIG. 1: Phase diagram in the plane of Temperature and Pinning Field in a schematic mean field system. At small temperatures 
for large values of the pinning field the system at equilibrium stays near the reference configuration (high overlap phase) while 
at small values of the pinning field the system stays far from the reference configuration (low overlap phase). The upper curve is 
line of dynamical transition below which a dynamically stable high overlap phase exist. The middle curve is the thermodynamic 
first order transition line between the low overlap and the high overlap phases. The lower curve is the spinodal of the high 
overlap phase, i.e. the point where it becomes unstable. The three curves terminate in the critical terminal point. (Figure 
from: S. Franz and G.Parisi, PRL 79 2486 (1997)). 



The effect of the frozen particles is similar to the one of an adsorbing matrix in a porous medium or of a pinning field, 
with the important addition that in this case the unfrozen particles remain in the original equilibrium state. 

The detailed phase diagram for spin models on random graphs was computed in |13( I14j . A theoretical discussion 
of the physical relevance of this situation for liquids and glasses and an exact computation for the mean field p-spin 
models were presented in |15) . Differently from the case of the pinning field where the field transforms the glass 
transition into first order and spinodal transitions, in the case of frozen particles the dynamic and thermodynamic 
transition lines keep with their glassy random first order character that one finds at zero pinning. 

In all these cases, the existence of a critical terminal point is interesting because while the dynamical critical line 
has to disappear in finite dimension thanks to dynamical activation, the critical terminal point, which is also the 
terminating point of the thermodynamic transition line could survive in finite dimensions and can be studied in 
numerical simulations and experiments. 

As it is usual for lines of phase transitions terminating in a critical point, the critical terminal point lies in a different 
universality class of the line. The critical properties of discontinuous dynamical transitions has been recently analyzed 
in JJ . It has been proposed that the time independent part of the fluctuations in the (3 and early a dynamical regimes 
admit a description in terms of a cubic replica field theory. The leading singularities of this theory in perturbation 
theory happen to coincide with the ones of a (p'^ field theory in a random magnetic field. At the critical point the 
coefficient of the (f>^ term vanishes and it is natural to make the hypothesis that the next relevant term is a (f>* term 
so that the resulting theory is the standard 0'*-RFIM {16j. 

Arguments in this direction have been put forward in [15] using a RG procedure. Unfortunately, the arguments 
in [IS], though suggestive, are not fully convincing. They are based on a Migdal-Kadanoff renormalization scheme, 
which uses a hybrid formalism where replicas are used to average out the randomness in couplings, but additional 
quenched disorder introduced to mimics the effect of the frozen particle is kept unaveraged. 

In this paper we use the tools of replica field theory to support the hypothesis that glassy critical points are in the 
universality class of the ^"'-RFIM. We analyze in detail the case of a dynamical transition line terminating in a critical 



point. Such a scenario applies exactly in the case in which a finite fraction of particles are pinned in an equilibrium 
condition. The case of a pinning field or of an adsorbing matrix presents additional complications that will be left 
to future work. Replica field theory can be used to find out the nature of this transition. We have to consider a 
system with n clone in presence of some constraint: when the number n of components goes to 1, one of the clones 
is at equilibrium and the other clones feel the effect of a quenched field. Since one does not specify which of the 
clones is privileged the final theory is replica symmetric. From a field theoretical perspective, at first site the RFIM 
hypothesis is self-evident: indeed if we take care of the leading terms, the replica theory corresponds to theory with 
a random temperature, that maps on a theory with a potential V{(f)) with a random magnetic field and the critical 
terminal point is described by a cj)'^ interaction. However this argument holds only for the leading terms and neglects 
sub-leading orders that may play a crucial role if the leading terms cancels. 

More precisely the mapping of the dynamical transition to the ^■^-RFIM comes from the extraction of the most 
singular contribution of a cubic replica field theory, with multiple fields of different scaling dimension. The neglected 
sub-leading terms turn out to become dimensionally relevant below dimension 6, the same dimension as the quartic 
terms of the ^''-RFIM. The contribution of cubic vertexes cannot therefore be directly dismissed on the basis of 
dimensional analysis. To understand the nature of the terminal critical point and the possibility of cancellations 
that restore the RFIM mapping, a careful analysis of the perturbative series is needed. The scope of this note is to 
investigate this problem at the level of Ginzburg criterion, computing the one loop corrections to the propagators due 
to the residual cubic vertexes and comparing them to the ones coming from the quartic terms. We find that a-priori 
unexpected cancellations are present so that the cubic vertexes contributions appear to be irrelevant as compared to 
the quartic ones. 

In order to test our theoretical results, we consider a diluted p-spin model (random XOR-SAT problem) on a random 
graph in presence of frozen spins and perform extensive numerical simulations to study finite size scaling close to the 
critical point. The simulations on large systems fully confirm our analysis, making us confident that the one-loop 
result indeed holds to all orders. 

In the next section we present the theoretical analysis. In the subsequent one the numerical simulation. We then 
conclude the paper. An appendix is devoted to the technical details of the theoretical computations. 

II. ANALYSIS 

Our starting point is the replica field theory ^7\ [18] 

-C = - / dx I ^(Vc/lab)^ -f TOi ^ 0^b + rn2 X! '/'a60ac + ^g ^ ^ab^crf 1 " g'^lTr 0^ - gt^2 ^1 '/'afa > (1) 

\ ab ah abc abed / ab 

where the field <j)ab{x) is an n x n space dependent matrix order parameter with vanishing diagonal terms, describing 
the fluctuations of the correlation function around its plateau value. Among all possible cubic replica invariant we 
have retained only the ones giving rise to the leading and possible next to leading singular behavior. It has been 
show in [5j that in the limit n — > 1 the theory describe fluctuations close to the dynamical transition, that occurs for 
mi — > (while m2 and rng remain finite). 

According to the discussion in [5] the components of the matrix (pab do not share a unique scaling dimension. 
Similarly to what found in the RFIM problem [TO], in order to use dimensional analysis one needs to change basis 
and define linear combinations (p, uj and Xab of the (pab that exhibit good scaling properties. The actual linear 
transformation reeds 

(Pab^ {(!)- ^^){^- Sab) + Uabl^ + Xab (2) 

where Uab — ^Q.a-(-i)» is a size two band matrix and Xab is a symmetric matrix null on the diagonal and such that 

^ Xa& = for all a 

b 

Y. UabXab = 0. (3) 

ab 

The analysis of the resulting quadratic form (see below) shows that 0, uj and Xab have well defined scaling dimension 
D^, D^ and D^ that, fixing the dimension Dmi = Ij verify D^^ — D^ + 1 and D^ = D^ + ^ Fixing the dimension 
of (j) by the condition that the action is adimensional one finds D^ — ^ ~ 1. Notice that if the normalization of 




FIG. 2: One-loop self-energy diagrams of a generic cubic theory. 

the dimension would be fixed by the condition that the dimension of the momentum is one, we should multiply all 
dimensions by a factor two. 

In terms of the new fields, keeping only the terms giving rise to the leading singularities in perturbation theory and 
sending n — > 1, one finds that (fTl) is equivalent to the Parisi-Sourlas action pO] 

dx -(to2 + m3)uj^ + u {~A(I) + mi(j) + igcf^) + ^'^Xab (-A + toi + 6.g0) Xab , (4) 

a.h 

where the matrix of fields Xab{x) has (n — f)(7i — 2)/2 — 2 independent components. These fields enter quadratically 
in (Hi). Explicit integration over them gives rise to a determinant to the power that tends to —2 for n -^ 1 and is 
equivalent to the one generated by fermion fields. The cubic vertex has a coupling constant g — lji—ijJ2, it has scaling 
dimension if — 2 and becomes relevant below dimension 8. 

In different systems the various parameters appearing in the action are function of the physical control parameters 
e.g. temperature and pinning field amplitude, or fraction of blocked particles. At the critical point one has wi — ui2^ 
the cubic coupling constant g is equal to zero and sub-leading singular vertexes need to be considered in the Landau 
free-energy expansion. We should therefore consider quartic terms and sub-leading cubic ones. Dimensional analysis, 
through repeated use of (pi), shows that among all possible quartic replica symmetric invariants, the most singular 
ones are [25] 

UiTr (^"^ + -U2 ^ (f)\i, + "3 X! ^lb<l>ac(l)cb (5) 

ab abc 



that after the change of basis the leading terms (lowest scaling dimension) give rise to the RFIM terms 



(m1 + U2 - U3) 



ab 



(6) 



The scaling dimension of this vertex is y — 3 which thus becomes relevant below 6 dimension. In addition one needs 
in principle to consider the cubic vertexes that do not vanish for uji = uj2- Among these, the ones of higher scaling 
dimension are Tr y^ and 'Y^^.b Xab; which for wi — L02 give rise to the term loi (Tr x'^ + Tlab Xah) • This term does not 
have evident reasons to vanish. Its superficial scaling dimension (when multiplied by dP x) is ^ — | and in absence 
of additional cancellations is would also become relevant below 6 dimensions. The relevance of these vertexes would 
spoil the equivalence with the RFIM at the critical point. 

In order to investigate the possible relevance of these terms in the action we could in principle analyze their 
contribution to the correlation functions of the various fields in perturbation theory. To do this one needs to keep 
into account the constraints (p| which give rise to some cumbersome combinatorics. 

We have preferred therefore to turn back to the original replica field theory (fTl) and compute the exact contribution of 
the cubic vertexes to the self-energy and to the correlation functions to the one loop level. The diagrams contributing 
to the self-energy of a generic cubic theory are depicted in fig. 2, they are the diagrams irreducible against cuts that 
separate the ending points PJ. 

The structure of the bare propagator of the theory in momentum space induced by replica symmetry is 



Gab:cdik) = (1 - Sab){l - 5cd) 



gi{k) 



Sac^bd + ^adSbc 



.92 (fc)- 



35d - 



.93 (fc) 



(7) 



where the three coefficients gi{k), i — 1,2,3 can be expressed as function of the mass parameters rrii appearing in the 
action ([l|, in the hmit ?i — > 1: 

9iik) = -^^ (8) 
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(mi + A;2)(-2(mi + fc^) + ^2) 
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(toi + fc2)2(2(TOi + k^) — 7712) 



(9) 
(10) 



Notice that while gi and 32 have a simple pole behavior for small k and mi, 173 has a double pole. 

Given the structure of the propagator and the vertexes, the direct manual evaluation of the self-energy diagrams 
to include the sub-leading order is rather awkward. The formal expression of the diagrams of fig. 2 is given in 
the appendix. Fortunately, the computation of the diagrams, which basically involves sums of constants and delta 
functions over replica indexes and momentum integration, can be fully automatized. We have achieved this goal using 
the algebraic manipulator Mathematica, defining functions that implement the Kronecker delta on the one hand and 
the sum of constants and delta functions on the other. In order to compute the contribution to the correlation function 
one needs to further multiply the resulting expression to the right and to the left by the bare propagator, which can 
also be done by automatic means. 

The result for the corrections to the correlation function can be cast in a form similar to the one of the bare 
propagator ([7]), parametrized by three numbers ri, r2 r^ analogous to the g's of ([7|. The resulting expressions, 
available upon request, are rather lengthy and we do not reproduce them here. For uji ^ uj2 the leading singularity 
behaves as m^ , with corrective terms of order m^ and higher. Comparing the leading contribution with 

the behavior of the bare propagator m^^ one sees that, as already found with the (/)^-RFIM mapping of pj, non 
Gaussian fluctuations become important below dimension 8. For wi = uj2 the leading terms vanish. Remarkably the 
computation shows that at the same time the corrective terms of order m^ , also vanish. In fact these contributions 

are proportional to cji — W2- The first non vanishing contributions are of order m-i . These are irrelevant down to 

dimension 4, as opposed to the quartic RFIM contribution which become relevant in dimension 6. 

This result proves the 0'*-RFIM mapping to the one loop order in perturbation theory. Though we believe the result 
to be valid to all orders the general mechanism for cancellation of the contribution of the cubic vertexes remains to 
be found. It would be tempting to think that for n —>■ 1 the fields Xab{x) in ([2]) verify fermionic algebraical relations 
that go beyond the simple identification of the replica determinant with the fermionic one that was used in |2j . 



III. SIMULATIONS 



In order to check whether the universality class of the terminal critical point is that of a ip'^-KFIM we have simulated 
the 3-spin model on a random regular graph with fixed degree z — 8. The Hamiltonian of the model is 

zN/3 



^ = - 2^ JtiSii^SqSii, 



where St — ±1 are N Ising variables and the zN indexes {ii,*2'*3} ^'^^ randomly chosen such that each variable 
appears exactly z times and each interaction involves 3 different variables. The couplings are independent random 
binary variables, generated according to P(J;^ = 1) = r and P(Jp = — 1) = 1 — r. In the thermodynamical limit, 
the model has been solved with the cavity method [22_ and presents a dynamical phase transition at a temperature 
Td = 1.3420(5) [2, followed by Kauzmann ideal glass transition at a lower temperature Tk- As long as T > Tk, 
the thermodynamical and dynamical properties of the model do not depend on the choice of the coupling bias r, 
because almost any coupling configuration has the same free energy, thanks to the fact that the annealed average 
equals the quenched one (logEjZj = Ej log Zj). In particular we can choose the couplings according to the Nishimori 
prescription [23] that corresponds to r = (1 -I- tanh/3)/2. This choice has the great advantage that the all-spin-up 
configuration Si = 1 is an equilibrium configuration and we can use it as the starting configuration for the study of 
the equilibrium dynamics and also as the initial configuration for a Monte Carlo simulation, which does not require 
thermalization .24.. 
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FIG. 3: The dynamical critical line as a function of the fraction 9 of frozen variables in the 3-spin model on a random regular 
graph with a fixed degree z = ?,. Critical temperatures are shown with error; the curve is just a guide for the eyes, not a fit. 



A. Phase diagram with frozen variables 



We are interested in studying the 3-spin model defined above, when a fraction 9 of variables are frozen to an 
equilibrium configuration (the all-spin- up configuration in the present case). The phase diagram in the {0,T) plane 
resembles closely the one derived in [T31II1] for the 3-XORSAT model 26J in the {6, a) plane (see Fig. 3 in 14J). It is 
also very similar to the one derived in [T5] for the spherical 3-spin model. 

In Fig. [3] we show the numerical data for the dynamical critical temperatures as a function of 9 obtained by the 
cavity method. All the critical points shown (but the terminal point that we discuss in the following) as been obtained 
by looking at the overlap with the reference equilibrium configuration [Si = 1): the dynamical critical temperature 
corresponds to the negative jump in this overlap when the temperature is slowly increased. 

A much more careful discussion requires the determination of the terminal critical point, that wc estimate to be 
located in 9c — 0.1534(2) and Tc = 1.5277(5). Indeed at the terminal critical point the phase transition is no longer 
of first order and there is no jump in the overlap with the reference configuration to be exploited (even for 9 < 9c but 
close to the terminal point the jump is so small that is not useful to estimate the transition location). On the contrary, 
for 9 > 9c the model has no phase transition at all and one can at most look at the precursors of the continuous phase 
transition taking place at the terminal critical point. 

The problem of locating the terminal critical point with high accuracy is delicate and we have followed a method 
based on solving the cavity equations in population, while looking at the instability parameter A. More precisely, 
we run the Belief Propagation (BP) algorithm for determining the fixed point distribution of cavity magnetizations 
P{m) that solves the equation 



/z—l 



m~ J 



n,(l + tpm\ml) - n,(l ~ tpm\ml) 
n,(l + tpm{m^) + n,(l - tpm{m^) 



where tp = tanh(/?). Then we have computed the local stability of the fixed point distribution by adding a small 
perturbation to each element m of the population and checking whether such a perturbation grows or decreases under 
BP iterations: the stability parameter A is defined such that the perturbation goes like e^* for large times t. So the 
fixed point is stable only if A < 0. In a discontinuous transition the condition A = identifies the spinodal points 
where the two states becomes locally unstable; while in a continuous transition the condition A = marks the unique 
critical point. 

In Fig.|4]we show the stability parameter A as a function of the overlap q with the equilibrium reference configuration. 
For each value of 9, the plot contains few tens of thousands of points measured at different temperatures around the 
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FIG. 4: The stability parameter A as a function of the overlap q with the reference configuration for two values of the fraction 
of frozen variables 6 = 0.512 and 9 — 0.156. Data points are from fixed point distributions measured at several different 
temperatures around the critical one. Curves are the best fitting quartic polynomials. Data at 6 = 0.152 show a discontinuous 
transition, while data at 9 = 0.156 show the absence of any critical point. In this plot, the terminal critical point 9c would 
correspond to a maximum A value equal to 0. 

critical one (both above and below the critical temperature). The noise in the data is due to the stochastic nature 
of the BP algorithm and to fluctuations related to the finiteness of the population (we have used populations of sizes 
10^). It is worth stressing that A mainly depends on g, while its dependence on the temperature is rather weak: 
indeed Fig. |4] shows data measured at different temperatures that lay on the same curve. The interpolating curves 
shown in Fig. [4] are quartic polynomials. 

For 9 = 0.152 the model undergoes a discontinuous phase transition: the two q values for which A = bracket 
an interval where no state can exist in the thermodynamical limit (since it would be locally unstable); by varying 
the temperature, the thermodynamic overlap has a jump not smaller than the size of this interval. Moreover the 
larger overlap where A = is the plateau value at the dynamical critical temperature T^. For 9 = 0.156 the stability 
parameter A never becomes positive and consequently the model has a unique paramagnetic phase and no phase 
transition. The 9c value of the terminal critical point is clearly in between these two 9 values and for 9 — 9c the X{q) 
function must have a maximum of height A = 0. The quartic interpolations shown in Fig. Hallow us to obtain reliable 
estimate of position and height of the maximum of the function X{q) for each 9 value. Fitting these maxima we then 
arrive at the following estimates for the terminal critical point: 9c = 0.1534(2) and qc = 0.499(1). 

A reliable way to estimate the temperature of the terminal critical point is by plotting the average overlap as 
a function of the temperature for 9 values at and around 9c (see Fig. pi). In this way we can estimate Tc as the 
temperature at which the average overlap measured at 9c reaches the critical value qc- The uncertainty on this value 
comes both form the statistical error in this fitting procedure and from the propagation of the uncertainty on the 
value 9c- In Fig. [5] we plot data for several 9 values such that one can appreciate how much the average overlap 
changes by varying 9. The final estimate for the temperature of the terminal critical point is Tc = 1.5277(5). The 
error is mainly given by the uncertainty on the value of 9c- 



B. Fluctuations at the terminal critical point 



Once we have obtained a reliable estimate for the location of the terminal critical point, 9c — 0.1534(2) and 
Tc = 1.5277(5), we have run extensive Monte Carlo simulations for those critical parameters. We have simulated sizes 
up to iV = 3 • lO"*. For each size, but the largest, we have simulated Iff* different samples (for the largest size only 
2500 samples were used). Thermalization is not an issue, given that the all-spin-up configuration is an equilibrium 
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FIG. 5: Average value for the overlap with the equilibrium reference configuration as a function of the temperature for several 
9 values. The lowest dataset shows a clear hysteresis cycle at the first order phase transition. The data plotted with a black 
line have been measured at the critical 9c- The horizontal line is the critical value for the overlap. 



configuration by construction. 

As already discussed in detail in P], there are 3 different sources of randomness in this model: the coupling 
configuration, the starting equilibrium spin configuration and the thermal noise. Thanks to the fact the annealed 
approximation is exact in this model above Tr-, in the thcrmodynamical limit every sample behaves exactly the same, 
when the average over starting configurations and thermal noise is performed. However, we do not take the average 
over many initial spin configurations, and we only use S'i = 1 as the starting configuration to avoid thermalization. 
So, the average over the many coupling configurations we consider does actually correspond to the average over the 
initial spin configurations: indeed we could gauge transform each sample in order to have roughly the same couplings 
and the initial spin configuration would change from sample to sample. 

So, we are left we only 2 sources of randomness: heterogeneities in the initial configuration (het) and thermal noise 
(th). Using the angular brackets for the thermal average and the square brackets for the average over the initial spin 
configurations, two different susceptibilities can be defined as follows 



Xth = N [{q') - {qf] , 
X,et = N{[{qf]^[{q)f) 



(11) 
(12) 



where the former measures thermal fluctuations within the same sample (averaged over the samples) , while the latter 
quantifies sample to sample fluctuations in the mean value of q (which is as usual the overlap with the equilibrium 
reference configuration). The total susceptibiUty is the sum of the two: Xtot = N {[{q'^}] — [{q}]'^) = Xth + Xhet- 



A dimensional analysis of action of the 
susceptibilities 



-RFIM leads straightforwardly to the following scaling relations for the 



Xth = N'^'xth{N'^'{q~qc)) , 
Xhet = N^^^Xhet^N'^'^iq-qc)) , 



(13) 

(14) 



In Fig. |6]we show the data collapse using the theoretically predicted exponents. We can see that the scaling for the 
largest fluctuations, Xhet, is very well verified, while some corrections to the scaling are still present in the data for 
the thermal fluctuations at the sizes we have simulated. 
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FIG. 6: Fluctuations over the thermal noise xth (above) and over the initial heterogeneities Xhet (below) rescaled according 
to the exponents predicted by a ^''-RFIM theory, for comparison in the inset we plot the same quantity using the scaling 
appropriate for the (/>^-RFIM. The data have been measured at the terminal critical point of the diluted 3-spin model by Monte 
Carlo simulations. 



IV. SUMMARY 



Summarizing, wc have found that Replica Field Theory of glasses predicts that fluctuations close to glassy critical 
points can be mapped into the Random Field Ising Model. This observation calls for the search of glassy critical 
points in liquid systems. In simulations one can use the method of pinning field or blocking a fraction of the particles. 
In experiments critical points could be found in liquids confined by porous media. 
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V. APPENDIX 

This appendix is devoted to the explanation of some technical details of the evaluation of the one loop corrections 
to the propagator. 

The diagrams that we need to evaluate are conceptually simple, but practically complicated by the presence of the 
replica indexes in the propagator and the vertexes that we write as: 



Gab-cdik) = (1 - (5a6)(l - Scd) 



9i(k) 7, +92{k) +93{k) 



(15) 



^ab-cd-ef = ^[^iSacSaeSbd^be + ^2SbcSdeSfa\ (16) 

where S denotes symmetrization against exchange of indexes in each couple (ab) (cd) and (e/) and exchange of the 
couples among themselves. The first self-energy diagram reads 

^abldik) = Jdq S (c.2^^^G<,e;c/(g)G,e;d/(fc + <z)+ (17) 

W1W2V' [Gae:cd{q)Gbe:cd{k + q) + Gab;ce{q)Gab:de{k + q)]+ (18) 

^ — ^e 

^lGab-cb{q)Gab:cb{k + q)) (19) 

(20) 

the second one is 

Sg.,(fc) = S/d<7 {^lT.eJ,gGef-jMGea-,acmM+U:lY.ejGe},eMGef-,abmac5bd + 

Wia;2 Y.eJs Gef:fg{q)Geg;ab{0)SacSbd + ^1^2 J2ef Gef-ef {q)Gef-ac{0)Sbd] ■ (21) 

In order to compute the correction to the propagator the one loop self-energy 'Slj^.^^[k) should be multiplied to the 
right and to the left by the bare propagator. 

SGab-Mk)= Y. Gab;a'b'{k)Ki'b';c'd'ik)G,,d'-cd{k)- (22) 

a' ,b' .c' .d' 

The whole calculation involves sums of constants and delta functions over replica indexes. In order to avoid mistakes 
we have automatized the calculation through the use of the software Mathematica. We reproduce here the commented 
worksheet used in the calculation. 
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Mathematica worksheet 
ln[i]:= (* Definition of Kronecker Delta Kd[a-b] *) 
ln[2]:= Unprotect [Times] ; Vnprotect [Power] 
Out[2]= {Power} 

ln[3]:= Kd[0] := 1 

JCd[-a_] .-= Kd[a] 
Kd[a_] Kd[a_] := Kd[a] 
Kd[a_] '■ 1 := Kd[a] 
Kd[a_] ^n_ := Kd[a] /; n > 1 
Kd[x_] := /; X f 

in[9]:= (* End Definition of Kronecker Delta Kd [a-ib] *) 

ln[10]:= 

(* Definition of sum *) 

in[ii]:= sum [Kd [a_ - b_] x_ , a_] := x / . a -* b 
sum [Kd [a_ - b_] x_ , b_] := x / . b -> a 
sum [Kd[a_ - b_], b_] := 1 
sum [Kd [a_ - b_] , a_] := 1 
sum [Kd[a_ - 1], a_] := 1 
sum [Kd[a_ - 2], a_] := 1 
sum [Kd [- 1 + a_ ] , a_ ] := 1 
sum [Kd [- 2 + a_ ] , a_] := 1 

sum[y_ x_ , a_] := y suni[x, a] /; FreeQ [y , a] 
sum [y_ + x_ , a_] := sum [x , a] + sum [y , a] 
sum [y_ , a_] := y n /; FreeQ [y , a] 

ln[22]:= (* End Definition of sum *) 

in[23]:= (* Propagator : to be considered a matrix in the couples of indexes {ab) and (cd) *) 

ln[24]:= G [a_ , b_ , c_ , d_ ] ; = 

Expand [{1 -Kd [a-i3]) (1 - Kd[c - d]) (gl / 2 (Kd[a-c] Kd[b - d] + Kd[a - d] Kd[b - c]) + 
g2 / 4 (Kd[a - c] + Kd[b - d] * Kd[a - d] + Kd[b - c]) * g3)]; 

ln[25]:= 

(* Increase the recursion limit *) 

ln[26]:= $RecursionLimit = 7000 

Cut[26]= 7000 

in[27]:= (* Matrix with the same structure of the propagator and different parameters *) 
ln[28]:= M [x , y , c , d] = G [x , y , c , d] I . {gl -^ ml, g2 -* m2 , g3 -* m3}; 

(* Matrix product *) 
in[29]:= Uno[a_, b_ , c_ , d_] := Evaluate [sum [sura [Expand [G [a, b, x, y] M[x, y, c, d]], x], y]] 

(* Notice that the product has the same structure of G and M *) 
ln[30]:= (* Find gl , g2 , gS such that G is the inverse matrix of M *) 

ln[3i]:= inv = FullSimplify [ 

Solve [{Uno[l, 2, 1, 2] == 1, Uno[l, 2, 1, 3] == 0, Uno[l, 2, 3, i] == 0}, {gl , g2 , g3)] / . n -* 1] 

r r 2 4 m2 2 {- 2 ml m3 + m2 (m2 + m3) ) ^^ 

cxjt[3i]= \\gl ^ — , g2 ^ , g3 ^ 

^ ml ml{-2ml + m2) ml^ (2 ml - m2) '' 

in[32]:= (* Introduce the momentum dependence ; k is adimensional *) 
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ln[33]:= invk = inv / . ml -* ml (1 + k '' 2) 

r r 2 4 m2 

Qjt[33]= gl ^ , g2 -) 

^ ^ [l + k^) ml [l + k^) ml [-2 (l + k^) ml + m2) 

2 [-2 [l + k^) ml m3+ m2 {m2 + m3)) 

g3^ 1} 



II + k^]^ ml^ l2 ll + k^] ml - m2] 



ln[34]:= 



ln[35]:= (* Representation in term of matrix 

elements : The replica symmetrix matrices have elements of three types : 
1) the couples of indexes coincide (ab) = (cd) 2) there is one 
index in common e.g. a = c but b ^ d 3) the indexes are all different. Calling xl, 
x2 , x3 these elements one can express gl,g2,g3 as a function of the x's *) 

ln[36]:= ss = Solve [{xl , x2 , x3}=={G[l, 2, 1, 2], G[l, 2, 1, 3], G[l, 2, 3, 4]}, {gl , g2 , g3}][[l]] 

at[36]= {gl ^ 2 {xl - 2 x2 + x3) , g2 ^ 4 {x2 - x3) , g3 ^ x3} 

ln[37]:= (* Same for the m's *) 

ln[38]:= ssl = SS / . {gl -* ml , g2 -* m2 , g3 -* m3 , xl -* yl , x2 -* y2 , x3 -* y3} 

CXjt[38]= {ml -^2 [yl - 2 y2 + y3) , m2 ^ 4 {y2 - y3) , m3 -^ y3} 

ln[39]:= (* ... and the elements of the matrix product, for n -* 1 *) 

uu = Evaluate [Simplify [{Uno[l , 2, 1, 2], Uno[l, 2, 1, 3], Uno[l, 2, 3, 4]} /. ss /. ssl I . n ^ 1\\ 

Qjt[39]= {2 [xl yl - 2 x2 y2 + x3 y3) , 2 [xl y2 - 2 x3 y2 + x2 [yl - y2 - 2 y3) + 3 x3 y3) , 
2 [4 x2 y2 + xl y3 - 6 x2 y3 + x3 [yl - 6 y2 + 6 y3) ) } 

in[40]:= (* Elements of product : expression above written explicitly *) 

in[4i]:= prod[xl_, x2_, x3_, yl_, y2_ , y3_] : = 

{2 (xl yl - 2 x2 y2 * x3 y3) , 2 (xl y2 - 2 x3 y2 + x2 (yl - y2 - 2 y3) * 3 x3 y3) , 
2 (4 x2 y2 + xl y3 - 6 x2 y3 + x3 (yl - 6 y2 + 6 y3))} 

in[42]:= (* End Matrix Product *) 



ln[43J:= 



(* Computation of diagrams *) 



ln[44]:= (* First Self -energy diagram *) 

(* Building up the wl - wl sub-diagram : Hll[a,b , c , d] *) 

in[45]:= h = Expand[G[a, c, al , cl] (G[b, c, bl , cl] I . {gl -♦ ggl , g2 -» gg2 , g3 -» gg3})'\; 
(* two internal propagators attached to wl - wl vertexes ; ggl, gg2, 
gg3 will be useful later to take into account momentum dependence *) 

in[46]:= Gil la_ , b__ , al_ , bl_] := Evaluate [(1 - Kd [a - b]) (1 - Kd [al - bl]) sum [sum [h, c], cl]] 
(* Sum over internal indexes *) 

in[47]:= Hll[a__, b_ , al_ , bl_] : = 

(Gll[a, b, al, bl] + Gil [b, a, al , bl] + Gll[a, b, bl, al] + Gll[b, a, bl, al]) I 
4 (* Symmetrization *) 

(* Building the w2 - w2 sub-diagram : H22[a,b,c,d] *) 

in[48]:= G22 [a_ , b_ , al_, bl„ ] ;= G[a, b, al, bl] (G[a, b, al, bl] I . {gl ■* ggl, g2 -* gg2 , g3 ■* gg3)) 
(* two internal propagators attached to w2 - w2 vertexes *) 

in[49]:= H22 [a_ , b_ , al_ , bl_] := G22[a, b, al , bl] 

(* Building the wl - w2 sub-diagram H12[a,b,c,d] *) 

in[50]:= p= sum [Expand [G [a, c, al , bl] (G[b, c, al, bl] I. {gl -* ggl, g2 ^ gg2, g3 ^ 993})], c]; 

(* two internal propagators attached to wl - w2 vertexes + sum of one internal index *) 
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in[5i]:= G12 [a_, b_, al_ , bl_ ] ;= Evaluate [p {1 - Kd[a - b])] 

in[52]:= H12 [a_ , b_ , al_ , bl_] := 1/ 2 (G12[a, b, al, bl] + G12[b, a, al, bl]) (* Symmetrization *) 

in[53]:= Diag[a_, b_ , c_ , d_] := wl ^2 Hll[a, b, c, d] + 2 wl w2 H12[a, b, c, d] + 

w2 ^ 2 H22[a, b, c, d] (* Total self -energy diagram (sura of the three) *) 

in[54]:= invkl = invk /. {gl -» ggl , g2 -* gg2, g3 -* gg3, k -* q} 

r r 2 4 m2 

ait[54]= \\ggl ^ , gg2 



ml {1 + q^) ml {l + q^j [m2 - 2 ml [l + q^)) 

2 [m2 {m2 + m3) - 2 ml m3 (l + g^)) 

gg3^ }} 

ral^ [l + q^)^ [-m2 + 2 ml [l + q^)) 

in[55]:= (* Substitution of parameters of propagator with inverse of raomentum dependent mass *) 

ln[56]:= H = Join [ invk [[!]], invkl [ [ I ] ] ] 

r 2 4 m2 
Cut[56]= gl -^ , g2 -^ , 

^ [l + k^) ml (-^ + *^) "-^ (-2 (^ + *^) ™^ + "^1 

2 (-2 (l + it ^) ml ra5+ m2 {m2 + m3)) 2 
g3 ^ , ggl ^ , 

(I + k')' ml' (2 [1 + Jt^) ml -m2) mi (i + q') 

4 in2 2 (ni2 (m2 + m3) - 2 ml m3 (l + g^) ) 
gg2 ^ , gg3 ^ 1 



ral (I + q') [m2- 2 ml (1 + <j^)) ^i^ h + ^2)^ (_m2 + 2 ml (l + g 
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in[57]:= Diagl [a_ , b_ , c_ , d_ ] = Diag[a, ib, c, d]/. ii; 

in[58]:= (* Multiply external propagators at zero raomentum to the self- energy *) 

in[59]:= (* Correction to the propagator : An integration over q is needed, 

remembering that q is a-dimensional and the integral should be multiplied by ml^{D/2) *) 

ln[60]:= F = Simplify [prod [G[l, 2, 1, 2], G[l, 2, 1, 3], G[l, 2, 3, 4], 

Diagl[l, 2, 1, 2], Diagl[l, 2, 1, 3], Diagl[l, 2, 3, 4]] I. (inv) /. n ^ 1][[I]]; 

ln[6i]:= dG = Simplify [ 

prod[G[l, 2, I, 2], G[l, 2, 1, 3], G\l, 2, 3, 4], F[[I]], F[[2]], F\\3\\\ I. (inv) I. 
n^l l.k ^-g][[l]]; 

ln[62]:= (* Series expansion of dG up to l/ml^4 terms *) 
ln[63]:= final = Series [dG , {ml, 0, -5}] 

16 {m2 + ra3)^ [5+4 q') {wl - w2)' 16 (ra2 + m3) (2 + q') {wl - w2) (5 wl - w2) 

Cut[63]= 



(l + g^)%l« (l + g^)'ml^ Olml]" 

16 (m2 + ra3)^ (5+4 g^) {wl - w2)' 16 (ra2 + m3) {wl - w2) il2 wl + 5 q' wl - 4 w2 - q' w2] 



{l + q')'ml' {l + q')^ml' 

1 16 (m2 + m3)' [5 + 4 q') {wl - w2)' 



0[ml]' (l + g2)%I« 

16 (m2 + m3) {wl - w2) il3 wl + 5 q' wl - 5 w2 - q' w2] 



[l + q'Yml' Oiml]'^ 

in[64]:= (* At the Critical point w2 = 

wl both the terms of order ml ^(-6) and the terras ml ^(-5) vanish *) 
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ln[65]:= final_critical = Simplify [Series [ 

Simplify [prod [G[l, 2, 1, 2], G[l, 2, 1, 3], G[l, 2, 3, 4], F[[l]], F[[2]], F[[3]]] I. inv / . 
n -, 1 /. W2-, wl /. k ^ q-[[[lU, {ml, 0, -4}]] 

12 wl^ 1 12 wl^ 1 12 wl^ 

Cut[65]= 



{l + q')'ml' 0[ml]' [l + q^)' ml' 0[ml]' [l + q^)' ml' 0[ml]'' 

in[66]:= (* The final result is that both the term l/ml^6 and l/ml^5 vanish for wl = 

w2 therefore the cubic vertexes can be neglected with respect to the quartic ones . 
Thee moment integration provides a factor ml^(D/2), 
so that the term is of order ml^(D/2-4), 
to be compared to the bare contribution ml ^(-2), 

telling that the terms start to contribute in D=4. The quartic terms on the other 
hand are of the order ml '^(D/ 2-5) which become relevant in dimension 6. *) 

(* The same is true for the second self- energy diagram *) 
in[67]:= (* Second Diagram : Ball *) 

ln[68]:= 

in[69]:= (* Ball without external legs *) 

in[70]:= ball [a_ , b_] = 

Expand [{ sum [G [c , a, c, b] , c] wl + G[a, b, a, b] w2) /. {gl -» ggl , g2 ^ gg2, g3 ^ 993}] 
(* First Ring *) 

gg2 wl gg2 n wl ggl w2 gg2 w2 1 

Cut[70]= 2 gg3 wl + + gg3 n wl + + + gg3 w2 - — ggl wl Kd [a - b] + 

2 4 2 2 2 

111 

gg3 wl Kd[a - b] + - ggl n wl Kd [a - b] + - gg2 n wl Kd [a - b] ggl w2 Kd [a - b] - 

2 4 2 

1 11 

- gg2 w2 Kd [a - b] - gg3 w2 Kd [a - b] + - gg2 w2 Kd[-a + b] gg2 w2 Kd [a - b] Kd [-a + b] 

2 4 4 

in[7i]:= (* Attach a propagator to the ball *) 

in[72]:= K [a_ , b_] = Evaluate [sum [sum [Expand [ball [c , d] G[a, b, c, d]], c], d]]; 

ln[73]:= 

(* Attach an external propagator with a wl vertex *) 

in[74]:= J [a_ , c_ , al_ , bl_] = sum [Expand [K [a , b] G[c, b, al, bl]], b]; 

in[75]:= (* Addition of external legs *) 

ln[76]:= Fl = wl Simplify [prod [G [1 , 2, 1, 2], G[l, 2, 1, 3], 

G[l, 2, 3, 4\, J[l, 2, 1, 2], J[l, 2, 1, 3\, J [1 , 2, 3, 4]]]; 

in[77]:= (* Attach an external propagator to the ball with a w2 vertex *) 

ln[78]:= L [a_ , b_ , c_ , d_'\ = Expand [K [a , b] G[a, b, c, d]]; 

in[79]:= F2 = w2 Simplify [prod [G [1 , 2, 1, 2], G[l, 2, 1, 3], 

G[l, 2, 3, 4\, L[l, 2, 1, 2], L[l, 2, 1, 3], L[l, 2, 3, 4]]]; 

in[80]:= (* Total ball diagram *) 

ln[8i]:= Ftot = Simplify [Simplify [Fl + F2] /. invkl I . inv /. n -^ 1]; 

(* This is the propagator correction up to a q integration *) 
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ln[82]:= Fser = Simplify [Series [Ftot , {ml, 0, - 4}] ] [ [I ] ] [ [I ] ] 

32 {ml + m3)^ (wl - v2)^ 8 (m2 + m3) [vl - w2) [4 m3 wl + m2 {{l3 + 6 q^) wl - [3+ 2 q^) w2) 

CXjt[82]= 



(I + q')^ ml' m2 [1 + q')^ ml' 

4 lie m3^ wl (wl - w2) + 

m2^ (1 + q^Y ml' 

4 m2 m3 i(20 + 15 q^ + 6 q') wl^ - [22 + 21 q^ + 10 q'') wl w2 + 4 [l + q^Y w2^\ + 

m2^ ((73+ 69 q^ + 24 q') wl^ -2 (39+ 45 q^ + 20 q') wl w2 + (l7 + 33 q^ + 16 q') w2^)\ + 

1 32 {m2 + m3)^ [wl - w2)^ 8 [m2 + m3) (4 m3 wl + m2 (5+2 q^) {3 wl - w2)) {wl - w2) 



0{ml\^ (1 + q^Yml' m2 (l + q')' ml' 

4 1 16 m3^ wl {wl - w2) + 

m2^ (1 + q^Y ml' 

4 m2 m3 i(22 + 15 q^ + 6 q') wl^ - (24 + 21 q^ + 10 q') wl w2 + 4 (l + q^Y "^^l + '^^ 



(3 (29+ 25 q^ + 8 q") wV -2 (47 + 49 q^ + 20 q") wl w2 + (19+ 35 q^ + 16 q") w2^) + , 

' [ml ] ^ 

32 {m2 + m3)^ {wl - w2)^ IS {m2 + m3) {wl - w2) (2 m3 wl + ml ((8+3 q^) wl - (3+ q^) w2)) 



(1 + q')^ ml' m2 (1 + q')^ ml' 

8 Is m3^ wl {wl - w2) + 

m2' (1 + q'Y "^^ 

2 m2 m3 {(23 + 15 q^ + 6 q') wl^ - (25 + 21 q^ + 10 q') wl w2 + 4 (l + q^)^ w2^\ + 

\ i 1 
m2^ ((47 + 39 q^ + 12 q") wl^ - (51 + 51 q^ + 20 q') wl w2 + 2 (5 + 9 q^ + 4 q') w2^) + 

' 0[ml]^' 
ln[83]:= Simplify [Fser / . w2 -* wl] 

16 {3 m2 + 2 m3) wl^ 1 16 {3 m2 + 2 m3) wl^ 1 16 {3 m2 + 2 m3) wl^ 1 , 



at[83]= 

m2 (1 + q^) ml" 0[ml]^ m2 (1 + q^) ml' 0[ml]^ m2(l + q^)ml' 0[ml]- 

in[84]:= (* We See that also for this term the result 

is 0(1/ ml ^4) and therefore irrelevant in dimension 6 *) 
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